home *** CD-ROM | disk | FTP | other *** search
/ Monster Media 1996 #15 / Monster Media Number 15 (Monster Media)(July 1996).ISO / prog_c / cuj0696.zip / DWYER.ZIP / LIB / IGAMI.C < prev    next >
C/C++ Source or Header  |  1996-03-30  |  2KB  |  105 lines
  1. /*                            igami()
  2.  *
  3.  *      Inverse of complemented imcomplete gamma integral
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * double a, x, y, igami();
  10.  *
  11.  * x = igami( a, y );
  12.  *
  13.  *
  14.  *
  15.  * DESCRIPTION:
  16.  *
  17.  * Given y, the function finds x such that
  18.  *
  19.  *  igamc( a, x ) = y.
  20.  *
  21.  * Starting with the approximate value
  22.  *
  23.  *         3
  24.  *  x = a t
  25.  *
  26.  *  where
  27.  *
  28.  *  t = 1 - d - ndtri(y) sqrt(d)
  29.  *
  30.  * and
  31.  *
  32.  *  d = 1/9a,
  33.  *
  34.  * the routine performs up to 10 Newton iterations to find the
  35.  * root of igamc(a,x) - y = 0.
  36.  *
  37.  *
  38.  * ACCURACY:
  39.  *
  40.  * Tested for a ranging from 0.5 to 30 and x from 0 to 0.5.
  41.  *
  42.  *                      Relative error:
  43.  * arithmetic   domain     # trials      peak         rms
  44.  *    DEC       0,0.5         3400       8.8e-16     1.3e-16
  45.  *    IEEE      0,0.5        10000       1.1e-14     1.0e-15
  46.  *
  47.  */
  48.  
  49. /*
  50. Cephes Math Library Release 2.0:  April, 1987
  51. Copyright 1984, 1987 by Stephen L. Moshier
  52. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  53. */
  54.  
  55. #include "mconf.h"
  56. #define    mtherr(a,b)
  57. extern double MACHEP, MAXNUM, MAXLOG, MINLOG;
  58.  
  59. double igami( a, y0 )
  60. double a, y0;
  61. {
  62. double d, y, x0, lgm;
  63. int i;
  64. double igamc();
  65. double ndtri(), exp(), fabs(), log(), sqrt(), lgam();
  66.  
  67.  
  68. /* approximation to inverse function */
  69. d = 1.0/(9.0*a);
  70. y = ( 1.0 - d - ndtri(y0) * sqrt(d) );
  71. x0 = a * y * y * y;
  72.  
  73. lgm = lgam(a);
  74.  
  75. for( i=0; i<10; i++ )
  76.     {
  77.     if( x0 <= 0.0 )
  78.         {
  79.         mtherr( "igami", UNDERFLOW );
  80.         return(0.0);
  81.         }
  82.     y = igamc(a,x0);
  83. /* compute the derivative of the function at this point */
  84.     d = (a - 1.0) * log(x0) - x0 - lgm;
  85.     if( d < -MAXLOG )
  86.         {
  87.         mtherr( "igami", UNDERFLOW );
  88.         goto done;
  89.         }
  90.     d = -exp(d);
  91. /* compute the step to the next approximation of x */
  92.     if( d == 0.0 )
  93.         goto done;
  94.     d = (y - y0)/d;
  95.     x0 = x0 - d;
  96.     if( i < 3 )
  97.         continue;
  98.     if( fabs(d/x0) < 2.0 * MACHEP )
  99.         goto done;
  100.     }
  101.  
  102. done:
  103. return( x0 );
  104. }
  105.